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This paper presents a new optimization procedure for Stirling engines based on neural network concepts. 
Based on modeling and experimental data an intelligent and fast method is proposed which finds the 
best values for different design variables. Design variables of Stirling engine are optimized using Multi- 
Layer Perceptron (MLP) neural networks. The optimization procedure is performed for three typical 
design variables for a given precision which has the capability to be extended for various types of engines 
and designs. 

Design variables are assumed to be phase angle , displacer stroke and working frequency of ST500 Stirling 
engine. Output variables which are subject to be optimized are the engine output power and efficiency. 
Maximizing a weighted average of these output variables is defined as total procedure objective. 

Usually, output variables are calculated using mathematical modeling based on structure, drive 
mechanism, dynamics and thermodynamic properties of the Stirling engine. The Multi-Layer Perceptron 
network is trained via 125 samples of design and output variables, generated by Nlog thermodynamic 
analysis code , in order to learn the relations between them. Then, the maximum point for the weighted 
average of the network outputs is determined with adequate precision. The proposed optimal point for 
the Stirling engine finally is substituted into the analysis code, and its outcome is compared with 
network results at this point. The proposed values are validated by analysis code. In order to check the 
validation of Nlog analysis code results for this scope of study, a test setup for ST500 engine is also 
established. Absolute pressure of working fluid and crank angle of the engine are measured instantly for 
three different charge pressures and compared with analysis results. 

© 2014 Elsevier Ltd. All rights reserved. 


1. Introduction 

Stirling engines have found a variety of applications as energy 
converters for Stirling/Dish, residential CHP, coolers and heat 
pumps etc. Specially, as renewable energy and energy efficiency 
have become more critical areas of research in recent years, the 
design and development of Stirling engine systems has been 
signified for combined heat and power systems. 


* Corresponding author. 

E-mail addresses: mazdak_hooshang@live.com (M. Hooshang), r.askari@ut.ac.ir 
(R. Askari Moghadam), s.alizadehnia@ip-co.com (S. Alizadeh Nia), m.t.masouleh@ 
ut.ac.ir (M.T. Masouleh). 

http://dx.doi.Org/10.1016/j.renene.2014.09.012 

0960-1481/© 2014 Elsevier Ltd. All rights reserved. 


The ideal Stirling cycle, as illustrated in Fig. 1, is an idealized 
thermodynamic cycle, involves four thermodynamic processes, 
including two isothermal and two constant volume processes [1 . 

The theoretical efficiency of a Stirling cycle operating between 
high temperature T e and low temperature T c is the same as that of a 
Carnot cycle. Heat is delivered to the cycle at T e , part of the heat is 
converted to work and the rest is removed as rejected heat at T c . 

Carnot revealed in 1820s that the maximum efficiency available 
from an engine working under this cycle depends only on the 
temperature difference in the cycle according to Eq. (1) which is 
called Carnot cycle efficiency. 

U = l-£ (1) 

I e 
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Nomenclature 

/ 

engine frequency (Hz) 



V 

engine efficiency 

U 

network input 

P 

engine output power (W) 

Y 

network output 

Or 

total heat rate rejected from the engine (W) 

? 

network output-layer weight 

T c 

heat rejection temperature (°I<) 

w 

network mid-layer weight 

T e 

heat absorption temperature (°K) 

W 

optimization weighting vector 

V f 

gas velocity(m/s) 

r 

engine performance 

P 

gas density(kg/m 3 ) 

i 

number of network mid-layer neurons 

n 

control-volume normal surface vector 

i 

neuron index (i = 1,2,...,/) 

X 

control-volume length(m) 

j 

number of network inputs 


wall area of control-volume (m 2 ) 

j 

network input index (j = 1,2,...J) 

S 

connecting surfaces between adjacent control- 

H 

number of network outputs 


volumes (m 2 ) 

h 

network output index (h = 1,2,...,//) 

1/ 

volume of control-volume (m 3 ) 

a 

sum of network input signals (m) 

T 

gas temperature inside control-volume ( K) 

T 

network learning rate 

P 

gas pressure inside control-volume (Pa) 

5 

network deviation factor 

T f 

gas shear stress from control-volume wall (Pa) 

Q 

total network training samples 

u 

gas specific internal energy (J/kg) 

z 

number of grid points along each input axis 

Qc 

net absorbed heat of control-volume gas (J) 

0 

Crank angle (Radian) 

Qe 

net injected heat to control-volume (J) 

(p 

piston-displacer phase angle (Radian) 

m w 

mass of control-volume wall (kg) 

u 

piston rod length (m) 

c w 

specific heat of control-volume wall (kj/kg) 

h 

displacer rod length (m) 

T w 

temperature of control-volume wall (°K) 

Cl\ 

piston gudgeon pin & surface distance (m) 

r 

working gas constant (kJ/kg.K) 

C*2 

displacer gudgeon pin & surface distance (m) 

c v 

gas specific heat for constant volume (kJ/kg.K) 

Q3 

displacer height (m) 

Re 

Reynolds number of control-volume gas 


piston gudgeon pin & crank axis distance (m) 

Re^[/ 

Reynolds number of control-volume gas 

^2 

displacer gudgeon pin & crank axis distance (m) 

P 

viscosity of control-volume gas (Pa.s) 

Cl 

piston cylinder head & crank axis distance (m) 

Nu 

Nusselt number 

C2 

displacer cylinder head & crank axis distance (m) 

k 

control volume index 

d 

displacer stroke (cm) 

I< 

working gas conductivity (w/m.K) 

r c 

Crank radius (m) 

t 

time (s) 

Af 

gas flow area (m 2 ) 





Fig. 1 . Ideal Stirling engine cycle processes. 


The efficiency of a real Stirling engine is also determined 
through measuring power output and rejected heat according to 
Eq. (2) as follows: [1] 


More than twenty sets of parameters like pressure, tempera¬ 
tures, engine working frequency, power piston and displacer ge¬ 
ometries, heat-exchanger geometries have direct effect on engine 
output power and efficiency and they complicate the design 
procedure. 


In other words, one of the important issues about designing a 
Stirling engine is to find an algorithm to find proper values of the 
engine parameters in order to optimize the engine performance [2]. 

In order to find the best design point, the simplest method 
consists in assuming all possible situations for engine specifications 
and checking the engine performance at each point. The main 
problem in this approach is the large number of unknown speci¬ 
fications. Thus calculating the effect of each parameter by means of 
convenient methods, needs a significant time which makes this 
method irrational. 

In principle, designing an engine starts with finding a set of 
parameters based on experiences and the next step is optimization 
for other parameters to improve the performance [3 . 

Toropov has utilized the simplified approximations of the 
original objective/constraint functions, obtained by the weighted 
least-squares method [4]. It used the information about function 
values gained at several previous design points in order to improve 
the quality of simplified approximations to the end of reducing the 
total number of calls for complex engineering systems analyses. 

Hirata has studied the effect of some heat losses and irrevers¬ 
ibility on the engine performance indices [5 . He has calculated the 
optimal performance and design parameters for maximum power 
and efficiency. It was considered that a compact Stirling engine was 
optimized by detailed measurements for heat balance. 

Timoumi et al. have carried out an optimization of GPU-3 engine 
data [6]. It resulted in reduction of the losses and to a notable 
improvement in the engine performance. The parameters of this 
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prototype were first applied on the developed model; the results 
were very close to the experimental data. Then, the influence of 
each geometrical and physical parameter on the engine perfor¬ 
mance and the exchanged energy of the regenerator were studied. 

Mahkamov et al. have developed a mathematical model of the 
LTD Stirling engine which accounts for hydraulic and heat losses in 
the working process and mechanical losses in the engine [7]. This 
model was optimized by Genetic Algorithm. As a result, a set of 
design parameters for the engine was obtained which provides a 
considerable improvement in the performance. Results were cali¬ 
brated using 3D CFD modeling technique. 

Chen et al. have introduced finite-time thermodynamics to 
optimize the performance of a solar powered Stirling engine [8]. 
The genetic algorithm is used to reveal the maximum output power 
of a Stirling engine by determining the thermal efficiency, the 
collector temperature and their corresponding values of the opti¬ 
mized parameters. The results of the simulation herein described 
can be directly applied to optimize the design of a practical solar 
powered Stirling engine. 

All of the discussed optimization algorithms are classified in 
gradient or evolutionary category. The gradient methods have a 
local view in the design variables space. The functional gradient on 
the point, defines the direction to get closer to optimum point. 
Depending on the starting point, after a few iterations, most of 
design variable space zones will lose their opportunity to contain 
the optimum point. Although these methods are fast, the final 
result dependent on the starting point commonly occurs whenever 
they are used for a nonlinear system, like a Stirling engine. It usually 
expressed as sticking to local optimums. Evolutionary algorithms 
such as genetic algorithm which use natural selection, mutation or 
recombination techniques to find the optimal point, reduce 
significantly the lack of sticking to local optimums but, in the other 
hand, they need numerous calls of calculation of engine perfor¬ 
mance which requires fast computer machines with significant 
time. 

Neural networks can be used for optimization purposes [9 . In 
Refs. 10,11 , neural network is applied to the end of predicting the 
power and performance of Stirling engines. For the purposes of this 
paper, this technique is used in order to define design parameters 
for optimizing the power and efficiency of Stirling engine. In fact, 
neural networks learn nonlinear and multivariable relations using 
large and rich set of data. Its main advantage consists in fast esti¬ 
mation of the output regarding to analytical algorithms and 
simulation codes. However, it should be trained for a new engine to 
have a suitable estimation of the output. 

In this article, a new optimization procedure based on neural 
network concepts is presented. It has been used for ST500 Stirling 
engine, designed and fabricated by IPCO. The network could predict 
the engine output power and efficiency in a fraction of a second 
with less than 2% error. Neural networks have high capability to 
approximate multivariable and nonlinear functions and it would be 
able to predict engine outputs and efficiency for other Stirling de¬ 
signs similar to ST500. It is obvious that new engine data should be 
learned by network, for each new engine. 

This is the first Stirling engine optimization procedure by means 
of a global view in the design variables space without needing 
numerous calls of calculation of engine performance. 

As the First step, the engine is modeled thermodynamically to 
calculate the engine output power and the efficiency at some 
sample points in design variable space. According to these samples, 


1 Low Temperature Difference. 

2 Gamma-type Stirling engine manufactured by IP-CO (Irankhodro power train 
company ). 


the neural network will learn the relations between engine per¬ 
formance and variation of design variables. At the end, the network 
predicts the engine performance for the whole design space with 
specified resolution and finds the optimal point. Global view of this 
method against local view of some of previously used methods 
reduces the possibility of sticking to local optimum points. More¬ 
over, the thermodynamic analysis takes significant time to calculate 
the engine performance for each sample point with enough accu¬ 
racy. Using this method a minimum distance between the sample 
points have been kept in the design variables space to prevent 
excess calls of the thermodynamic analysis for areas of the space. 
This feature makes this method rather fast. 

Another characteristic of this method is fast estimation of a 
Stirling engine performance, which is a benefit of using the neural 
network. Once the network is trained, it could be used in order to 
estimate the engine performance for specified conditions in frac¬ 
tion of a second. This benefit comes in handy for control purposes 
of FPSE type which would be our next research step. It would help 
if the controller knows about the amount of energy being produced 
at any time. In Section 2, design variables are selected according to 
this purpose. 

To determine the thermodynamic behavior of Stirling engine 
working fluid, there are a number of Stirling cycle analysis tech¬ 
niques. The techniques are categorized from zero to fourth order 
analysis. A zero order analysis is an experimental relation between 
engine output power and some of design parameters. 

In a first order analysis, engine output power is calculated 
assuming ideal isothermal expansion and compression. For sinu¬ 
soidal volume variations a closed form solution is also available 
[ 12 ]. 

Second order analysis is an extension of first order in which 
some power losses due to limitation of heat transfer and gas 
pressure drop and finally subtracted from isothermal calculation. 
Recently, Cheng [13] proposed a second order analysis which is 
suitable for prediction of thermal efficiency of a beta-type Stirling 
engine. 

A third order analysis is a one-dimensional fluid modeling of 
working gas inside the engine to calculate power generation, heat 
transfer losses and pressure drop effects simultaneously. 

Future detailed models like 3D modeling of working gas are 
classified in fourth order Stirling engine analysis. Makhkamov and 
Djumanov have presented a full 3D calculation of a Stirling engine 
using their CFD model in Ref. [14 . 

In this paper, calculation of engine performance is done using 
Nlog thermodynamic analysis code which is a third order analysis. 
The code results verified by experimental data and its error is less 
than 18% for the cycle work. Then using the foregoing code, a set of 
data is produced. Accordingly, the neural network trained to 
recognize relations between inputs and outputs of the dataset. In 
the next step, the whole design variable space will be searched 
using the network to find the best optimum point. The procedure is 
done for ST500 Stirling engine, with specifications shown in Table 1. 

Fig. 2 shows the ST500 Stirling engine used for the experiment. 
Heater, regenerator, cooler and compression space can be seen in 
this figure. 


3 Free Piston Stirling Engines, type of Stirling engines in which piston and dis¬ 
placer are not mechanically connected to each other and freely oscillate inside their 
chambers. In this type, an external magnetic field releases the power generated by 
the engine and through this way controls the piston oscillation amplitude and 
phase on a specific frequency. 

4 Stirling engine third order analysis code, provided by Stirling machine 
researcher M. Hooshang( mazdak_hooshang@live.com) 
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Table 1 

ST500 engine characteristics and operating conditions. 


Specifications 

Values 

Type 

Gamma 

Electrical output power 

500 W 

Total efficiency 

8.5% 

Standard charge pressure 

8 bars 

Working fluid 

Helium 

Working frequency 

14 Hz 

Fuel 

Natural gas 

Cooling substance 

Water 

Power piston stroke 

75 mm 

Power piston diameter 

84 mm 

Displacer stroke 

75 mm 

Displacer diameter 

97 mm 

Phase angle 

90 0 

Heater type 

Tubular (6 mm dia.) x 20 

Cooler type 

Duct (13 mm 2 section area) x 144 

Regenerator material 

Steel matrix (0.96 porosity) 

Heat absorption temperature 

350-420 °C 

Heat rejection temperature 

30-50 °C 


engine working frequency. Increasing or decreasing engine fre¬ 
quency may also changes optimum phase angle [1 . 

Although there are many other variables affecting on the engine 
performance, these three variables are enough for a typical opti¬ 
mization procedure to converge to an optimal point while other 
engine specifications remain constant. 

Both of the optimization variables represent the engine perfor¬ 
mance but to make a compromise, a weighted combination of them 
is defined in Eq. (Variables rj and p are normalized values of 
efficiency and output power respectively). 


r = w 


v 




Consequently the optimization variable vector will be reduced 
to a scalar form, which is dependent to the weight vector of W. The 
weight vector is chosen arbitrary considering the fact that sum of 
the vector elements must be equal to unit. 


2. Design & optimization variables 


3. Design variable ranges and required precision 


The target of this optimization process is to find the best value 
for design variables to keep the optimization variables as high as 
possible. Some of the engine parameters are assumed to be design 
variables and some others to be optimization variables. Generally, 
the target would be maximizing two variables (H = 2 ), engine 
output power and efficiency. Thus they are defined as optimization 
variables. Depending on the defined target, three design variables 
(J = 3) are selected including displacer stroke (d), engine frequency 
(f) and phase angle (<p). 

Changing the displacer stroke affects on displacer swept volume 
and compression ratio. Compression ratio is the main design 
parameter for a Stirling engine to be adapted with wall temperature 
of engine heat exchangers. This parameter is also easy to apply on 
the real gamma-type engine against displacer bore or power piston 
bore and stroke. Changing the swept volume will change maximum 
fluid flow rate and velocity which is related to pumping loss in the 
heat exchangers. There are different ways to control the flow ve¬ 
locity but the easiest and effective way would be changing the 



Heater 

Regenerator 

Cooler 

Flywheel 


Compression space 


Fig. 2. Main parts of ST500 Gamma type Stirling engine (Design and manufactured by 
IP-CO). 


Required precision for determination of displacer stroke and 
phase angle is based on geometrical tolerances of engine manu¬ 
facture. Required precision for determination of engine frequency 
also depends on dynamometer accuracy. 

A reasonable range for all design variables are also determined 
from experience and are shown in fable 2 . 


4. Neural network 

The total optimization procedure consists of two parts; pattern 
recognition and interpolation. 

First, consider the engine model nature as a static transfer 
function. The engine model gets the design variables and gives 
output power and efficiency. While the inputs change, each of the 
outputs may uniformly ascend, uniformly descend or consecutively 
ascend and descend. There is also a situation in which some of input 
variations negate the others; thus the outputs may remain almost 
constant. For example, convective heat transfer coefficient in a heat 
exchanger pipe increases whenever the piston swept volume in¬ 
creases. But it may remain constant if engine working frequency 
decreases simultaneously. 

Here is the role of the neural network to form such transfer 
function which can cover any of four discussed situations or even 
combination of them. Thus, a kind of pattern recognition procedure 
must be used which can support all these situations. Neural net¬ 
works could act as a beneficial pattern recognition tool in which 
each neuron transfer function acts as a part of - discussed effects. 

Above discussion lead us to a ridge-shape network transfer 
function. Structure of a sigmoid function, one of the well-known 
ridge-shape transfer functions, is shown in Eq. (4). In the forgoing 
equation, 05 is weighted summation of input signals as formulated 
in Eq. (5) and j/j is output signal from neuron i ( Uj is network input 
and Wji is mid-layer weight). 


Table 2 

Design variable ranges and required precision. 


Design variable 

d: Displacer 
stroke (m) 

6: Phase 
angle (deg) 

/: Frequency 
(Hz) 

Minimum acceptable value 

0.02 

60 

12 

Maximum acceptable value 

0.12 

120 

22 

Required precision 

±0.001 

±0.25 

±0.25 
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u \ 

u 2 

u 3 



Fig. 3. MLP neural network model structure with one mid-layer, three inputs and two 
outputs. 


yi = (1 +exp(-<7;)) 1 

( 4 ) 

J 


°i = J2 w fi u i 

( 5 ) 

j'=i 



power and efficiency. For each sample point, it takes at least 20 min 
using a machine with a T4500 CPU 4 GB of RAM. 

According to trial and error, z = 5 seems to be adequate for all 
axes. (z\ = z 2 = z 3 = z). Thus the network is trained viaq = 5 3 = 125 
samples. The total time spent for Nlog code to calculate power and 
efficiency of the sample points was about 50 h using the afore¬ 
mentioned machine. 

The two inputs and three outputs of each sample are normalized 
according to the following expression to be fed into the network: 


X 


X-min(X) v , r 

- TT7Z - . /vV X = P, 77, a, <pj 

max(X) - min(X) 



6 . Neural network training 

A well-known method to train a MLP neural network is error 
back-propagation (BP) algorithm [16]. 

BP is an algorithm which decides to set the network weights (? 
and /jl) such that the network output error is kept as low as possible. 
Network output error is defined considering engine output power 
and efficiency for each training sample as shown in what follows: 

e,=p-Y, (9) 


In this paper, as presented in Fig. 3, a Multi-Layer Perceptron 
(MLP) neural network with one mid-layer and the sigmoid transfer 
functions for the mid-layer neurons is used. MLP is one of the best 
beneficial tools which is widely used for pattern recognition and 
classification 15]. 

According to Eq. (6), the network output is weighted summation 
of neuron output values. Therefore, while the network inputs in¬ 
crease (U\,U 2 and U 3 ), depending on signs of f^s for each of neu¬ 
rons, the network outputs (Yi and Y 2 ) may increase, decrease, 
consecutively increase and decrease or remain constant. These are 
the four situations which the engine may respond. 


y/. = Ef»M ( 6 ) 

i=l 

A full input-output formula of this network can be formulated 
as follows: 


i 


Vh = £ffc 
1=1 


M + exp 




The network outputs are prediction of the network for engine 
output power and efficiency according to design variables dis¬ 
cussed before. 


e\ =rj-Y 2 (10) 

in which Y\ and Y 2 are defined to be the network outputs. 

Output power and efficiency of the engine are calculated using 
Nlog code for each value of design variables. 

Network inputs are also normalized values of design variables 
where one has: 

U,=d,U 2 = cp,U 3 =f (11) 

At each step of training, network deviation factor is calculated 
according to BP algorithm using Eq. (12). Then required changes for 
the network weights are calculated using Eqs. (13) and (14). The 
learning rate in these two equations determines the speed of 
weights changes. Assuming a greater value makes the network 
weights to change more quickly (Quicker pattern learning) but may 
cause the weights to change unstably and prevent from a good 


convergence. 


0(Tl h= 1 

( 12 ) 

^ih = T e h Yi 

(13) 


5 . Sample generation 

In Section 2, three design variables are considered (J = 3) which 
makes the design variable a 3-dimensional space with 3 axes. The 
optimization objective is to determine these variables such that the 
maximum of output power and efficiency is obtained. On each axis, 
a number of grid points (z) are defined with equal distances. 
Combination of grid points creates total 3-dimensional sample 
points which are used to train the neural network. The total 
number of the sample points is equal toq = 7 } = z 3 . The engine 
power and efficiency at each sample point is calculated by the Nlog 
code (as explained in Section 8 ). Increasing the number of grid 
points on each axis, will raises the accuracy of network prediction 
along that direction but there is time limitation in calculation of 


Aw^ = y <5; Uj (14) 

New values for the network weights are set using following 
equations: 

+ A? ih (15) 

wjf w - W" ld + A Wjf (16) 

The training procedure is done with 5 epochs. At each epoch the 
total of q = 125 training samples are applied to the network so the 
total steps of network training is 5 x 125 = 625. As discussed 


5 A 2.3 GHz Dual-Core type CPU manufactured by Intel Inc. 
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before, the network has three inputs (j = 1,2,3) and two outputs 
(h = 1 , 2 ). 

Mean squared error during network training epochs is shown in 
Fig. 4 assuming 10 neurons for the network mid-layer and y = 0.7 
for the network learning rate. 

In general, it is recommended to use simplest structure with 
minimum number of neurons to achieve fast interpolation with 
more reliability. Fig. 5 shows the variation of error norm for each 
network outputs according to the number of network mid-layer 
neurons. Thus 8 neurons for this layer makes the average error 
minimum and it is suitable for a better interpolation. 

7. Determination of optimal point 

To cover the full range of design variables space with precision 
specified in Table 2, estimation of engine output power and effi¬ 
ciency for 1.2 x 10 5 points is needed ( 0.12 - 0 . 02 / 

2 x 0.001 x 120 - 60/2 x 0.25 x 22 - 12/2 x 0.25 = 1.2 x 10 5 ). 
Therefore the total 1.2 x 10 5 points dataset is fed to the network 
and the network predicts corresponding outputs for each one. 
Engine performance (P) which is weighted combination of the 
network outputs is calculated for each point using weight vector of 
W = [0.7 0.3]. These weight values are assigned due to importance 



Fig. 5. Error norm according to number of mid-layer neurons. 


Table 3 

Neural network estimation comparison with Nlog code results for the achieved the 
optimal point. 


Optimization variables 

V 

P 

Neural network estimation 

13.21% 

878.2w 

Nlog code result 

13.45% 

880.6w 

Neural network estimation error 

1.78% 

0.27% 


of Stirling engine working efficiency against output power but any 
arbitrary vector with summation of elements equal to unit could be 
assigned instead. Next the global maximum of T is determined. 
Denormalization of the network outputs for the maximum point 
results estimated engine output power and efficiency for the 
optimal point which are shown in Table 3. Calculation of the engine 
output power and efficiency on optimal point using Nlog code is 
also shown in Table 3. 

Achieved optimal point for ST500 engine is. 

d opt = 10.2 cm, <p opt = 92.7 deg, / opt = 18.4 Hz 

with output power and efficiency of: 

Popt — 880.6 W, ?7 0 pt = 13.45% 

while according to Table 1, previous operating point of the engine 
was: 

dstd = 7.5 cm, <p std = 90.0 deg, / std = 14.0 Hz 
with output power and efficiency of: 

Pstd = 500 W > Pstd = 8 -50% 

Figs. 6-8 illustrate the variation of engine performance against 
two of three design variables while the third variable set to cor¬ 
responding optimal point value. 

Increasing 380 W of output power and 4.95% of efficiency shows 
the effectiveness of the optimization approach and low estimation 
error of neural network shows the considerable ability of MLP 
neural network to predict Stirling engine performance. 



Fig. 6. Neural network output for engine performance against engine frequency and 
phase angle (d = 10.2 cm). 
















































M. Hooshang et al. / Renewable Energy 74 (2015) 855-866 


861 


0 

o 

c 

0 

E 

o 

0 

Cl 

0 

c 

'o) 

c 

LU 


0.15 

0.1 

0.05 

0 


-0.05 

40 



I uu 


0.1 

Stroke(m) 


0.2 


Phase angle(deg) 


Fig. 7. Neural network output for engine performance against variation of stroke and 
phase angle (f = 18.4 Hz). 



and energy equations respectively [17 . But they are simplified to 
the one-dimensional form which are fit for this third-order 
analysis. 



(17) 


/ \ 



s k 


Free flow area for the working gas transformation is assumed to 
be constant for each control volume thus Eq. (17) is converted to Eq. 
(20). The area is also assumed to be constant during time which 
leads to obtain Eq. (21 ). 

/ PVf nAA f + l t (p k \x k ) = 0 (20) 

Sfc 


/ pvftidAf + A fk ^ (, Pk x k ) =0 (21) 

s t 


Density and velocity of the working gas are assumed to be 
uniform over cross sectional area between adjacent control vol¬ 
umes and also inside each control volume. Thus following equation 
becomes the final simplified continuity equation used in the code. 

E Ps v f s n sA fs + Af t ^ ( P k x k ) = 0 (22) 

s=S k 


Fig. 8. Neural network output for engine performance against variation of frequency 
and stroke (cp = 92.7 deg). 


8. Nlog code description 

In this work, learning patterns are generated via Nlog thermo¬ 
dynamic analysis code. Nlog is a third-order Stirling cycle analyzer 
which calculates the generated power and rejected heat by Stirling 
engine and also instant working gas state variables according to 
following parameters of the engine: 

1) Geometry properties of all gas transferring channels, pipes, and 
compression and expansion chambers. 

2) Geometries of linkages between engine moving parts. 

3) Engine initial charge pressure and local temperatures. 

4) Temperature of heat exchangers wall (Constant during time 
period). 

The code divides the whole gas transferring channels and pipes 
of the engine, to control volumes and determines gas dynamic and 
thermodynamic parameters for each control volume by solving 
continuity, momentum and energy equations simultaneously. 
Equations (17)—(19) are standard forms of continuity, momentum 


Extending above assumptions of working gas velocity, density 
and cross sectional area, for working gas pressure in Eq. (18) will 
result in Eq. (23). The following is the simplified form of the mo¬ 
mentum equation. 

0 

gf ( Pk A f k x k% ) + X P s n s A fs + E ( v fsPs) v f s n s A f s + T ft A w k = 0 

s=S k s=Sk 

(23) 

Energy equation is also simplified upon assuming gas specific 
heat to be equal for everywhere of each control volume and 
expressed as follows: 



The total heat transferred into each control volume can be 
calculated using the following relation: 


dQcfc _ Nu k KkA Wk , _ T \ 

d t ~ D h ( Wt k > 
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Variation of wall temperature for each control volume is 
calculated from: 

^ + —L_ (^ = 0 (26) 

dt c w m Wk \ dt dt J 

Gas specific internal energy is assumed to be only function of 
gas temperature according to Eq. (27) for Helium gas 18]: 

u k = c v T ]o (c v = 3.11 kJ/kg.K) (27) 

Gas thermal conductivity is assumed to be function of gas 
temperature and pressure according to Eq. (28) for Helium gas 18]. 

I< k = 2.6810 -3 x (l + 1.123P k 10- 8 ) (r^ 1-2 * 10 ’ 9 ^ ( 28 ) 

Thermodynamic Ideal gas behavior is considered and stated in 

Eq. (29) [18]. 

Pk = PiJTio 0r = 2.08kJ/kg.K) (29) 


Two assumptions are used by Nlog code depending on wall-type 
of each Stirling engine chambers. The assumption for a regenerative 
type walls is shown in Eq. (30-a) and for isothermal type is shown 
in Eq. (30-b). 


dQe t 

dt 



(30-a) 


—^ = 0 (30-b) 

dt 

Formulations of Reynolds number (Re) and kinetic Reynolds 
number (Re w ) are described in Eq. (31) to Eq. (34). 



(31) 





ft * 2 


t. 

A 




Fig. 9. ST500 linkage kinematics. 



8-k 


(32) 


^ = 1 - 475 x 10 "7tSo 


(33) 



(34) 


Fig. 9 represents the Kinematic modeling of the engine linkage. 
Referring to Fig. 9, <p, r c , l\, l 2 , d, ci, C 2 and a\ to a 3 are engine 
construction parameters, i.e., design parameters, and are constant. 
Following equations relate these parameters to b\ and £> 2 . 


l\ = r 2 c + b\ + 2 r c £q cos(tf) 


(35) 


/2 = y 2 c + b\ + 2 r c b 2 cos (6 + (p) 


(36) 


Upon solving above equations for b\ and £> 2 , they are expressed 
as functions of crank angle in Eq. (37) and Eq. (38). Thus x\ to X 3 
could be expressed as functions of crank angle according to Eqs. 

(39) and (41). 


b 


1 



1 

2 


r c cos 6 


(37) 



( r 2 cos( 2 q> + 28) 

V 2 



) 1 

2 

+ r c cos (<p + 0 ) 


(38) 


X] = C] — ci] — b ! ( 39 ) 

x 2 = c 2 -a 2 - b 2 (40) 

x 3 = d - a 3 - x 2 (41) 

For each control volume indexed by /<, Eqs. (22)—(29) contain 12 
unknown variables including x, r, Nu, v/, p, P, u , Qe, Qc, K T w and T. 
Depending on control volume index, the variable x is calculated 
from Eqs. 39 and 40 or (41) for variable volume chambers or set to 
constant for constant volume chambers. The values of r and Nu are 
expressed as functions of Reynolds number (Re) and kinetic Rey¬ 
nolds number (Re w ) according to experimental correlations have 
been presented by Kays and London 19], De Monte et al. [20] [21] 
and S. Y. Zheng [22 for different ranges of Re and Re w , type and 
dimensions of cross sectional area and wetted surface roughness 
for the engine heat exchangers. 

Eqs. (22) and (30) totally form a system of 9 equations which 
Nlog code solves to determine 9 remained unknown variables at 
each time step for each control volume. The whole engine is divided 
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Table 4 

Engine chambers and model control volumes. 


Engine chambers 

Control volumes 

Piston 

1 (k: 1) 

Connecting pipe 

6 (k: 2-7) 

Under displacer 

1 (k: 8) 

Cooler 

7 (k: 9-15) 

Regenerators (x5) 

6 (k: 16-45) 

Heater 

7 (k: 46-52) 

Above displacer 

1 (k: 53) 

Buffer 

1 (k: 54) 


to 11 chambers and 54 control volumes as shown in Table 4. A 
typical form of two adjacent control volumes is shown in Fig. 10. 
Runge Kutta 4th order method is used to solve the equations. 

Finally, the results including engine output power, p, and 
rejected heat, Q r are calculated as follows: 



Fig. 11. Test results for condition C. 


P= £ A fA- c w ( 42 > 

lcVariable 

volumes 




Fig. 12. Nlog code result for condition C. 
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Fig. 10. Thermodynamic variables of two adjacent control volumes. 


Fig. 13. Graphical comparison between test and Nlog code P-V diagrams for condition 
‘A’. 
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Fig. 14. Graphical comparison between test and Nlog code P-V diagrams for condition 
‘B’. 



Fig. 15. Graphical comparison between test and Nlog code P-V diagrams for condition 
‘C’. 


Qr= E de k ( 43 ) 

k: Cooler 
volumes 

The engine efficiency is calculated by Eq. (2) for the network 
training. 

9. Thermodynamic analysis results validation 

We need to make sure the Nlog code results to be reliable for this 
scope of study for optimization of ST500 engine. Thus a test setup is 
established for the engine that measures critical state variables for 
some random situations to validate the code results. 

For the validation, comparison for all the state variables of the 
test results and the code results along one time period could be 
done. But there are limitations for measuring some of the state 
variables like instant temperature of working gas which needs very 
low response-time sensors or piston position sensor which needs 
significant resistance against high temperatures for the sensor. 

Thus the piston chamber gas pressure and engine crank angle 
were assumed to be instant test results because it is possible to 


measure them and rather easy to read by general data acquisition 
units. Other instant variables like the temperature and chamber 
volumes could be expressed as functions of gas pressure and crank 
angle according to state equation of ideal working gas (Eq. (29)) and 
geometries of engine linkages (Eqs. (35) & (36)). 

A piezo-type pressure sensor is placed at the top of the engine 
compression chamber. These types of sensors are vulnerable to 
high temperature and this is where the gas temperature is less than 
any other chambers. A cooling adapter coupled to a water cooling 
system is also used to keep the sensor body temperature in 
acceptable range. 

A crank angle encoder is coupled to the engine flywheel. The 
encoder is a relative type but with trigger teeth it is used as an 
absolute type. According to tooth number, the crank angle and total 
working volume of the engine is calculated. 

Average value of the engine rejected heat also could be calcu¬ 
lated simply by measuring coolant liquid flow, coolant entrance 
temperature and coolant exit temperature. 

The whole test setup for ST500 engine is shown in Fig. 16. 

All of boundary conditions during the test were kept near to 
engine standard specifications presented in Table 1. 

The test is done on ST500 engine for three different charge 
pressures (Represented as initial conditions) and heat exchanger 
wall temperatures (Represented as boundary conditions) which are 
shown in able 5. 

P-V diagram for condition ‘C is shown in Fig. 11 which is a 
graphical result of the test. 

All of the ST500 engine specifications ( Table 1) are put into the 
Nlog code. Initial conditions of working gas pressure and boundary 
conditions of heat exchanger wall temperatures are set same as the 
test rig condition (Table 5). 

The code is run on a machine explained before. Total spent time 
for the machine was about 70 min 50 min to simulate engine warm 
up period (Equivalent to three complete cycles), and about 20 min 
to obtain main results (Equivalent to one complete cycle). P-V di¬ 
agram for condition ‘C is shown in Fig. 12 which is a graphical result 
of the code. 

Comparisons between the test and code results for three 
different conditions ( Table 5) are shown in Table 6. As shown in this 
table, the maximum estimation error is about 11 % for power and 
about 18% for rejected heat rate. All power losses of dynamometer 
system (Including the belt, magnetic field and electric losses of 
generator wiring) are pre-calculated and compensated by data 
accusation software. The code doesn't take into account mechanical 
losses due to bearings and rings. In order to have a rough estima¬ 
tion of these losses, friction loss value of the engine is extended for 
higher frequencies and added to experiments results values. It 
should be noted that due to the non-oil lubrication of the engine, it 
is possible to consider linear behavior for bearing friction factor. 

This prediction error level is low enough for an optimization 
procedure comparing with some other Stirling cycle analysis codes. 
For example, the errors of same comparison between uncali¬ 
brated version of GLIMPS code calculations and experimental 


6 Instant pressure sensor of Kistler 4075A. 

7 Compact cooling adapter of Kistler 7505B. 

8 Crank angle encoder of 365-series. 

9 Test which is performed by the engine manufacturer to determine internal dry 
friction losses of it at low frequencies. 

10 Only results of uncalibrated version of the codes is valuable for a design or 
optimization procedure. 

11 Stirling cycle analysis code being used by National Aeronautics and Space 
Administration (NASA). 
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Pressure sensor body 
Sensor coolant pipe 

Sensor data cable 

Encoder base 
Engine coolant pipe 



Flywheel belt 
Encoder connector 
Encoder 

Starter power supply 

Sensor temperature 
regulator unit 

Generator power cable 
Encoder data cable 

(a) 


Data accusation 
board 

Data accusation 
power 



State variables 
online monitoring 



Fig. 16. (a) ST500 engine test setup, (b) Data accusation unit. 
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measurements of RE-1000 engine are 1.5% for rejected heat and 
19.0% for power [23]. 

Comparisons between overall curve shapes of P-V diagrams for 
three different conditions (Table 5) are shown in Figs. 13 to 15. 

There is also an adequate compliance for overall curve shapes of 
P-V diagrams. Area inside the P-V diagrams represents the power 
exchanging between engine pistons and working gas. 

10. Conclusion 

In this article, an optimization procedure for Stirling engines 
using neural network was presented. Choosing MLP neural network 


A 1 kW free piston Stirling engine designed by Sunpower Inc. 


for this kind of pattern recognition problem was because of the 
ridge-shape transfer functions nature which enables the network 
to predict engine performance fast enough and precise. It was 
shown that via 8 mid-layer neurons, the network could predict the 
engine output power and efficiency in a fraction of a second with 


Table 5 

Three different initial and boundary conditions. 


Condition 

A 

B 

C 

Charge pressure (bar) 

8.1 

6.2 

4.6 

Coolant liquid flow (lit/min) 

6.1 

6.1 

6.1 

Coolant entrance temp. ('C) 

33.5 

40 

35 

Coolant exit temp. ('C) 

41.0 

46 

42 

Heater tip temp. ('C) 

393 

400 

395 

Heater base temp. ('C) 

353 

350 

355 


12 
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Table 6 

Comparison between Nlog code and experimental measurements of ST500 engine. 


Condition 


A 

B 

C 

Nlog Code 

Rejected heat rate (kW) 

3.331 

2.974 

2.270 


Power (kW) 

0.440 

0.363 

0.262 

Experiment 

Rejected heat rate (kW) 

2.94 

2.52 

2.52 


Power (kW) 

0.400 

0.343 

0.235 

Error 

Rejected heat rate (%) 

13.2 

17.8 

9.9 


Power (%) 

10.0 

5.8 

11.4 


less than 2% error (compared to the thermodynamic code result) 
while the thermodynamic code takes at least 20 min for one 
simulation session. Defining more design variables obviously takes 
into account more relations between inputs and outputs, so the 
number of the neurons must be increased to keep estimation error 
at low levels. Spacing of training samples along each input axes, 
also affect on network prediction quality; but according to trial and 
error it is determined that 5 points for each axis is enough to cover 
the whole common ranges of many design variables. 

Comparison between engine output power and rejected heat for 
test setup and Nlog thermodynamic analysis code in same condi¬ 
tions, represented acceptable compliance, so the code results could 
be accepted. Confirmation of 878 W of output power and 13.21% of 
efficiency by the code comparing 500 W and 8.5% for the engine 
before the optimization, illustrates the advantageous of optimiza¬ 
tion procedure using the neural network. 
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